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ABSTRACT 

The intergalactic magnetic field (IGMF) can be indirectly probed through its effect on electromagnetic 
cascades initiated by a source of TeV gamma-rays, such as active galactic nuclei (AGN). AGN that 
are sufficiently luminous at TeV energies, "extreme TeV blazars" can produce detectable levels of 
secondary radiation from Inverse Compton (IC) scattering of the electrons in the cascade, provided 
that the IGMF is not too large. We review recent work in the literature which utilizes this idea to 
derive constraints on the IGMF for three TeV-detected blazars-lES 0229+200, 1ES 1218+304, and 
RGB J0710+591, and we also investigate four other hard-spectrum TeV blazars in the same context. 
Through a recently developed detailed Monte Carlo code, incorporating all major effects of QED and 
cosmological expansion, we research effects of major uncertainties such as the spectral properties of 
the source, uncertainty in the UV - far IR extragalactic background light (EBL), undersampled Very 
High Energy (VHE; energy > 100 GeV) coverage, past history of gamma-ray emission, source vs. 
observer geometry, and jet AGN Doppler factor. The implications of these effects on the recently 
reported lower limits of the IGMF are thoroughly examined to conclude that presently available data 
are compatible with a zero IGMF hypothesis. 

Keywords: cosmology: large-scale structure of universe, cosmology: cosmic background radiation, 
gamma rays: general, magnetic fields, methods: numerical, radiation mechanisms: non- 
thermal 



1. INTRODUCTION 

Numerous observations have established the presence 
of magnetic fields in our own Galaxy and in the galax- 
ies, clusters, and filaments that comprise the large scale 
structure of the Universe. However, an unambiguous 
detection of the intergalactic magnetic field (IGMF), 
presumed to exist in the void regions, which repre- 
sent a significant fraction of the volume of the uni- 
verse, remains elusive. Such a field could be produced, 
for example, through astrophysical mechanisms such as 
bul k outflows of ma g netized material from r adio galax- 
ies (jKronberd (|1994f ): iKronberg et~aTl (|2001[ )). although 
it is unclear whether such proc esses cou l d effic iently fill 
the entire volume of the voids (Zweibel (2006)). Alter- 
nativel y, processes such as the Biermann battery mech- 
anism (jBiermannl (|195Qh ) operating during phase tran- 
sitions in the early universe could produce the IGMF, 
provided that its correlation length is sufficiently large to 
over come magnetic diffusion and survive to the present 
day (Gr asso fc Rubinstein! (j2001f )). "Primordial origin" 
hypotheses, such as this one, are particularly attractive 
because the IGMF could then play the role of the seed 
field necessary in magnetohydrodynamic models com- 
monly invoked to explain t h e fiel ds observed in galax- 
i es an d clusters (jWidrowl (|2002f ): iKulsrud fc Zweibell 
(|2QQ8f )). Consequently, the detection of the IGMF could 
provide important insights for solving outstanding prob- 



lems of its origin and role in both the cosmology and 
astrophysics of structure formation. 

The standard observational technique used to detect 
weak magnetic fields in galaxies, measuring the Faraday 
rotation of light from distant quasars, is inadequate for 
detecting the IGMF for two reasons. The first is that 
Faraday rotation measures the integrated magnetic field 
along the line of sight and therefore the determination of 
the IGMF relies on the subtraction of the imperfectly 
measu r ed Galactic magn etic field (jKronberg fc Perrvl 
(|1982[ ): iBlasi et al.l {1999)). The second, and perhaps 
more important one, is that a sufficiently weak IGMF 
strength will produce a Faraday rotation measure that 
is below the resolution limit of currently employed tech- 
niques. Existing Faraday rotation measurements place 
an upper limit of B < 10 -9 Gauss on the strength of 
an IGMF with a correlation length greater than 1 Mpc, 
and this limit weakens as the correlation length decreases 
until the limits due to Zeeman splitting measurements 
of absorption lines in dista nt quasars, become more cq n- 
straining (as summarized in lNeronov fc Semikoz (2009)). 
As a result of these two effects, until recently, only upper 
limits on the IGMF strength have been established. 

A new measurement technique has emerged during 
the past few years which may become a more sen- 
sitive tool for the measurement of IGMF characteris- 
tics. This technique relies on observations of blazars, 
believed to be AGN whose jet is oriented along the 
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line of sight to Earth, in the gamma-ray energy range 
from 100 MeV to greater than 10 TeV and is de - 
scribed by several authors ([Ne ronov & Semikoz (2006); 
lEung wanichavapant & A haronianl (j2009l ): iDolag et al.l 
([2009D : lEhdv et al.l ([2009D ). Briefly. TeV-scale gamma- 
rays from the blazar interact with the EBL, producing an 
electron-positron pair. The electrons and positrons then 
undergo IC scattering on the Cosmic Microwave Back- 
ground (CMB) radiation, producing secondary gamma- 
rays of a lower energy than the primary. As a result, an 
electromagnetic cascade develops. Because the pairs' tra- 
jectories depend on the Lorentz-force interactions with 
the magnetic field, the temporal and angular profiles of 
this cascade emission at the GeV scale carry information 
on the strength of the magnetic field in which the cas- 
cading occurs. If cascading develops in the void regions, 
then spectral, temporal, and angular characteristics of 
the secondary radiation will depend on IGMF proper- 
ties. A number of studies have characterized the spectral 
properties of the secondary photons as well as the tem- 
poral p rofile, commonly re fe rred to as an " echo" (jPlagal 
(fl995h : llchiki et al.l (|2008h : iMurase et all (120081 ) ), an d 
the angular profile, or "hal o," (lAh aronian et al.l ()1994f ): 
iNeronov fc Semikozl (j2009|) ; |Ah]irs| ([2011|)). Comparison 
of the characteristics of secondary gamma-ray radiation 
with existing data can therefore become the methodology 
for studying properties of the IGMF. 

Most published studies of the subject have focused 
on the spectral properties of the secondary radiation. 
For insta nce, utilizing sim p le geo metri c models for the 
cascade, Neronov Vovkl (|201Q[ ) and iTavecchio etTal] 
(2010) demonstrated that a lower limit on the IGMF 
strength can be placed by requiring that the secondary 
gamma-ray GeV-band emission does not exceed current 
measurements. In another study, ITavecchio et al.l (j2011l ) 
performed detailed modeling of the spectral energy dis- 
tribution of four blazars, thereby reducing the depen- 
dence of the conclusions on assumptions a bout the prop- 
erties of individual sources. Furthermore, iDermer et al.l 
(|2011l ) relaxed the assumption that the characteristic 
time to build up the secondary gamma radiation is less 
than the duration of the Very High Energy (VHE; en- 
ergy > 100 GeV) activity of the source, and derived a 
less constraining lower limit on the IGMF. Expanding 
on these simplified geometric models, Huan "eTatl (12011) 
described a semi-analytic model employing the energy- 
dependent distributions of electron and positron ener- 
gies in the cascade and accounted for the effects of the 
source lifetime. M onte C arlo simulat i ons w ere used by 
IDolag et all (|2011f ) and iTavlor et all (j2011[ ) to confirm 
and improve the results of si mplified geometric semi- 
analytic models. In particular, ITavlor et al.l (|2011[ ) em- 
ployed a three-dimensional particle tracking simulation 
to follow the cascade development and derived a lower 
bound on the strength of the IGMF which appeared to 
be consistent for three blazars studied. 

This paper reviews IGMF constraints derived from 
previous studies part i cularl y focusing on the conclu - 
sions of ITavlor -et al.l ([201 lh. [Neronov fc Vovkl fl2010l) . 
ITavecchio et al.l ([20111 ). and IDermer et al. | ([20lTh . Be- 
cause of the scientific importance of these results, we wish 
to systematically understand all the ways in which these 
IGMF constraints may fail. Particular emphasis is given 
to the effects of major uncertainties such as the spectral 



properties of the source, history of VHE activity, geomet- 
rical characteristics of the source, and uncertainty due to 
the poorly known EBL spectral energy density. The con- 
clusions are derived based on the newly developed Monte 
Carlo code which is described in section [2] together with 
models for the EBL, IGMF, and AGN source model. Sec- 
tion [3] describes the data utilized in this paper, while sec- 
tion [4] provides detailed analysis of individual sources and 
a comparison of the results of prior publications. Section 
O the discussion, concludes the paper with a brief review 
of alternative interpretations of the data. 

2. NUMERICAL SIMULATIONS 

VHE photons escaping a source such as an AGN jet 
interact with surrounding diffuse photon fields and gen- 
erate electromagetic cascades. Cascading occurs in mor- 
phologically complex environments of photon and mag- 
netic fields of the host galaxies, the large scale structure 
filaments, voids, etc. The amplitudes of the magnetic 
fields and the density of the photon fields in these struc- 
tures sensitively affect the temporal, angular, and spec- 
tral evolution of the cascading, secondary photons. 

For example, the highest energy of the escaping pho- 
ton is determined by the spectral energy density of back- 
ground photons of the host galaxy. A 30 TeV photon 
progagating through a Milky Way-like galaxy will have 
an optical depth of ~ 1, based on rough estimates of 
the energy density of the Galaxy in the far infrared 
(~ 100 jim). Photons with energies higher than this 
will either be absorbed by interactions with the galac- 
tic light or the CMB, on spatial scales less than the size 
of the galaxy (< 1 Mpc). These photons will initiate 
cascades under the influence of galac t ic m agnetic fields 
(~ 10- 5 - 10- 6 G, see, e.g. iWidrowl (120021) ) which are 
strong e nough to isotropize the s econdary photons of the 
cascade ([Aharonian et al.l (|1994f )). 

Photons with energies low enough to escape the galaxy 
are expected to predominantly interact with the EBL. 
If this interaction occurs in the environments of either 
galaxy cluste rs with m agnetic fields of order 10 -6 — 10 -7 
G (see, e.g. Wi drowl ([20021 )) or intervening large scale 
structure filaments, with magnetic fields > 10 -12 G, 
then the secondary electrons will be isotropized, thereby 
dramatically attenuating the observable flux of the sec- 
ondary photons produced by them. Effects of the electro- 
magentic cascading may become observable with present 
day instrumentation when the interaction with the EBL 
occurs in the voids with characteristic magnetic fields 
< 10- 12 G. 

In order to explore in detail the potentially observable 
spectral, angular, and temporal effects of cascading in 
the voids of the large scale structure, we have developed 
a fully 3-dimensional Monte Carlo code. It propagates 
individual particles of the cascade in a cosmologically 
expanding universe and accounts for all QED interactions 
with the EBL and CMB without simplifications. In the 
following, details regarding these numerical simulations 
are provided and some of its capabilities are illustrated 
with several characteristic examples. 

2.1. Cosmology and EBL Model 

Throughout this paper, a Friedmann Robertson 
Walker cosmology is used with critical matter, cosmolog- 
ical constant, and radiation densities given by Qm = 0.3, 
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Figure 1. Compa rison between the EBL model of 

Domingue z et all (2011) and the 3 EBL models used in the 
present work-a typical EBL model (EBL model 1), one chosen 
with a low dust peak (EBL model 2), and one chosen with a low 
stellar peak (EBL model 3). 

= 1.0 - Hm, and Q R = 8.4 x 10" 5 respectively. The 
radiation density is negligible for nearby sources of z < 1. 
The standard EBL model chosen in the simulations at z 
= is modeled with 6 energy density points (I\(240/jLm) 
= 15.9 nW m- 2 sr" 1 , J A (60) = 7.3, J A (12) = 2.0, 
Ia(2.5) = 9.4, 7 A (0.38) = 3.2, J A (0.12) = 0.43), and a 
cubic spline is used to find the EBL e nergy density a t 
intermediate wavelengths as adopted in IVassilievI ([2000). 
Figure [1] shows a comparison of this EBL model (labeled 
EBL model 1) with one recently proposed, which is based 
on various empirical observational d ata as well as on 
EBL m odels existing in the literature (jDommguez et all 
(|2011f )). Two additional EBL models are displayed, one 
(labeled EBL model 2) with a substantially lower dust 
peak, and another (labeled EBL model 3) with a reduced 
stellar peak. These models are explained in detail and 
used in section B~3l 

In this study, we investigate sources with redshift < 0.2 
and thus, we neglect evolution of the EBL in the co- 
moving reference frame; only the effect of cosmologi- 
cal expansion on the EBL energy density and energy of 
EBL photons is taken in to account. By comparison with 
iDommguez et al.l (|2011[ ), we do note that evolutionary 
effects of the EBL may need to be taken into account for 
sources of z > 0.3. 

2.2. Magnetic Field Model 

The observational effects of magnetic fields in the 
gamma-ray signal of extragalactic sources originate in 
the deflection of the charged particles from the trajectory 
of the primary photon and subsequent deflection of the 
direction of the secondary IC scattered photons. Differ- 
ent regimes of influence of the IGMF can be determined 
by exploring the interplay between the characteristic co- 
herence lengths of the IGMF, the e^ IC cooling length, 
and the distance from the interaction point to the ob- 
server. For the production of secondary photons above 
100 MeV, which is of relevance to this paper, the pair- 
produced e ± should have energies on the scale of > 100 
GeV. In the simulation, electrons are tracked down to en- 
ergies of 75 GeV. A 1 TeV electron loses its energy to IC 
scattering on the CMB over a characteristic length of 0.4 
Mpc, Xic oc E~ x . For the coherence length of magnetic 
fields, Xigmf ^> A/cs the IC scattering effectively hap- 
pens in the environment of nearly constant Bigmf- For 



Xigmf <^ A/c, the deflection of charged particles occurs 
in the non-coherent regime and leads to smaller deflec- 
tion. In this study, the coherence length of the IGMF is 
conservatively chosen to be 1 Mpc, which corresponds to 
coherent scattering for all photons of interest (with ener- 
gies > 100 MeV). It has been observed that the reversal 
field length of the magnetic fields in clusters of galax- 
les is on the scale of 10 - 100 kpc (|G rasso & Rubinstein 
(2001)) and reflects the spatial scales of the distribution 
of plasma. Thus, the magnetic fields in the voids with sig- 
nificantly larger characteristic plasma distribution scales, 
are likely to have coherence lengths much larger than 
this. 

The IGMF is modeled in the code as a system of cu- 
bic cells with a size equal to the coherence length of the 
IGMF and magnetic field amplitudes which are equal 
in value but randomly oriented in direction. To preserve 
cosmic variance, each cell is assigned an orientation when 
the first particle of the electromagnetic cascade propa- 
gates through it. If the cascade develops over a large 
number of these magnetic field cells, the observable ef- 
fects of the IGMF are randomized. However, if the dis- 
tance to the observer is comparable to the size of the 
magnetic field domain, the observational effects of the 
randomly chosen direction become significant. For this 
study, we analyze sources at distances greater than a few 
hundred Mpc (z > 0.1). The evolution of the IGMF is 
unknown, but is likely dominated by cosmological ex- 
pansion for z < 1. Therefore, the size of the domains 
and the magnetic field value are evolved with the stan- 
dard (1+z) -1 and (1+z) 2 dependences. The values of the 
magnetic field reported in this paper, refer to the values 
at z = 0. 

Figure [2] illustrates the mean time delay of secondary 
photons produced by a monoenergetic beam of 100 TeV 
primary photons at a redshift of z = 0.13. The pho- 
tons arriving at the observer are integrated over an aper- 
ture radius of 10.0°. For each of the energy bins (4 per 
decade), the mean delay time is computed for 6 magnetic 
fields including the zero field case. The time delay for the 
latter is due to QED scattering of the secondary particles 
of the pair production and IC scattering processes. For a 
10 GeV photon, the time delay amounts to about a half 
hour and for 0.1 GeV photons, the time delay is about 
10 hours. The saturation effect at low energies is due 
to the aperture cutoff. Figure [2^ shows the mean time 
delay with no EBL photo ns as IC targe t s in t he simula- 
tions. It is consistent with lTavlor et all (|2011[ ) Figure 2, 
and follows the spectrum of T delay E~%, derived ana- 
ly tically by making several si mplifications, as explained 
in iNeronov fc Sem ikoz (2009). Figure EJd shows the re- 
sult when the EBL field is included. The time delay for a 
non-zero field of secondary photons with energies above 
10 GeV is significantly increased, because electrons can 
move farther from the position of pair production and 
still scatter higher energy EBL photons towards the ob- 
server increasing the average time delay in a given en- 
ergy bin. The effect can be seen more clearly in Figure 
[3l which shows the distribution of arrival times of sec- 
ondary photons in a single energy bin of 10.0 - 17.78 
GeV with and without the target EBL photons for IC 
scattering. 
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Figure 2. Mean time delay of secondary gamma-rays produced by a monoenergetic beam of 100 TeV photons at a redshift of z = 0.13, 
for a varying Bjgmf with (a) (left) the EBL excluded as a target photon field in the IC scattering of electrons and (b) (right) with the 
EBL included as a target field. 



10 4 



- EBL included 
-EBL not included 




._ 10 3 10° 

rime Delay [yr] 

Figure 3. Distribution of arrival times of secondary photons in 
a single energy bin of 10.0 - 17.78 GeV at B IGM f = 10~ 17 G 
with and without the EBL photons included as a target for the IC 
scattering computation. 

2.3. Gamma Ray Source Model 

The gamma-ray source model employed in the sim- 
ulations is based on leading theoretical speculations 
about the nature of th e TeV blazar source (see, e.g., 
lUrrv fc PadovanH (|1995f )). In the reference frame of the 
blazar jet, the VHE photons are distributed isotropically 
with a power law spectrum of index a. Once this dis- 
tribution is boosted into the reference frame of the host 
galaxy with a Doppler factor of T = (1 — /3 2 ) -1 / 2 , it can 
be parameterized as 



e^f(e,0 v ) = F o S 2 (±- 



-a+l 



-e/e c 



where e is the photon energy in the reference frame of 
the host galaxy, S — [r(l - /3cos# v )] -1 , #v is the viewing 
angle from the blazar jet axis to the line of sight of the 
observer, F Q is a flux normalization factor, and e Q is an 
energy scale factor. To account for absorption of photons 
at the highest energies inside the host galaxy, an expo- 
nential cutoff at energy e c is introduced. This model is 
characterized by 4 independent physical parameters, and 
is sufficient to model the VHE part of the spectrum (> 



100 GeV). 

Since the energy range of interest for this study covers 
more than 5 decades of energy (from 0.1 GeV to >10 
TeV) the HE part of the source spectrum (< 100 GeV) 
is allowed to obey a power law with a different spectral 
index 7 



dF v ,2 j 

de 



(els) 



-7 + 1 



exp 



-a+l 



exp 



(e c £ ) cbS 
(=*) -S- 



< 1 



> 1 



(1) 



The sixth independent parameter, e#, introduced in 
this equation, is the spectral break energy. This multi- 
parameter gamma-ray source spectrum is given at the 
redshift of the host galaxy and is necessary and suffi- 
cient to satisfy observational data of TeV blazars in both 
the HE and VHE regimes. 

To illustrate the effects of different gamma-ray source 
model parameters, simulations are shown for z=0.13, 
B=10- 15 G, e c = 30 TeV, a = 7 = 1.5. Figure H shows 
the observed spectra of a source with the above parame- 
ters when it is viewed at different angles and different jet 
Doppler factors. All photons arriving within an aperture 
of 5° around the source are integrated. The black (solid) 
line shows the spectrum of the prompt photons reach- 
ing the observer, which is modified by absorption on the 
EBL. It is assumed that this observed spectrum is the 
same for different viewing angles and different Doppler 
factors. Figure H^i shows the spectrum of secondary pho- 
tons when the source is viewed at V = 0°, 2°, 5°, 10° 
and Figure HJd shows the spectrum of secondary photons 
for T = 5, 10, 30, 100. Viewing a source with the same 
prompt spectral energy density (SED), but with increas- 
ing viewing angle or Doppler factor implies the power in 
the jet must increase. 

The higher energy end of the secondary photon spec- 
trum (> 20 GeV) is unchanged for different viewing an- 
gles and jet Doppler factors because photons at these 
energies are generated from primary photons leaving the 
source with nearly zero deflection from the angle to the 
observer and the flux of these photons is directly pro- 
portional to that of the prompt photons. The lower en- 
ergy end of the spectrum (< 20 GeV) is generated by 



5 




10 ,_ 10" im 10" 
Energy [GeV] 



10 1<r 10 

Energy [GeV] 



1 

A 



2E 



de 



n e (e) 



Figure 4. Simulations of a source at z=0.13, with a = 7 = 1.5, e c = 30 TeV, and Bigmf=10 15 for a) (left) F 
viewing angles and b) (right) V = 5° and four different jet boost factors. 

secondary electrons of lower energies, the trajectories of 
which are significantly deflected from that of the pri- 
mary photon producing them. Additionally, the cooling 
length due to IC losses increases inversely proportionally 
to energy, allowing significantly larger deflection angles. 
The position of the peak of the SED is correlated to the 
value of the magnetic field, and its intensity is propor- 
tional to the overall power in the jet which increases with 
larger viewing angle or Doppler factors. When the char- 
acteristic angular size of the jet becomes larger than the 
viewing angle (4 > V ), the SED is nearly independent 
of the Doppler factor. 

Figure [5] displays a simulation of the photon arrival 
distribution from a source at z = 0.13, with gamma-ray 
source model parameters a = 7 = 1.5, e c = 30 TeV, 
B/gmf = 10 -15 gauss, T = 30, at four different observ- 
ing angles, V = 0°, 2° 5°, and 10°. The main trend 
in these figures is that the overall luminosity of prompt 
and secondary emission rapidly declines as the observing 
angle increases, and the the photon distribution around 
the source becomes increasingly axially non-symmetric, 
when 6 V ~ Y- The luminosity of the secondary photons 
relative to the prompt emission rapidly declines with in- 
creasing viewing angle, thus, detecting non-axially sym- 
metrical halos around AGN with existing instrumenta- 
tion may prove challenging. Figure [5] is in qualitative 
agreement with pr eviously report e d fin dings in the study 
of these effects bv lNeronov et al.l (j2010l ). 
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2.4. Particle Propagation and Interaction 

The development of cascades in intergalactic space 
intiated by the VHE photons is modeled by full 3- 
dimensional ray tracing and interaction of all particles: 
electrons and photons. The processes of pair production 
and IC scatter ing are treated with the use of fu l l QED 
cross sections (peli^l (HHH); IGould fc SchrideTI (jl967)) 
implemented without simplification. For example, both 
the EBL and CMB are included in IC scattering as seed 
photon fields and scattering includes the Klein Nishina 
regime. Both pair production and IC scattering are 
treated in the code through the sampling of marginal 
probability density functions. For example, the mean 
free path A of electrons of energy E due to IC scattering 
is given by 



(2) 

where /3 is the speed of the electron, cr ie (x) is the differ- 
ential IC cross section, x = (2Ee/m 2 e )(l — f3 cos 6) where 6 
is the collision angle and n e (e) is the spectral energy den- 
sity of the isotropic seed photon field. The code samples 
the propagation length of the particle with the use of the 
mean free path, A to determine the position of the inter- 
action. It then generates the marginal probability den- 
sity for an interaction of the electron with a given energy 
photon and samples the energy of the interacted photon. 
The process is continued by sampling the interaction an- 
gle of a given energy photon by generating a marginal 
probability density of the angular distribution of seed 
photons. The azimuthal angle is sampled randomly from 
a uniform distribution. At the completion of this process, 
the kinematics of the interaction is fully determined, al- 
lowing the computation of the energy-momentum vector 
of the outgoing photon and electron. The pair produc- 
tion process is treated similarly. 

Propagation of both photons and electrons is simulated 
in a cosmologically expanding universe. Energy losses of 
photons are only due to cosmological expansion, while 
energy losses of electrons also include IC cooling. All 
photons with energies above 100 MeV are tracked un- 
til they reach the z=0 surface, at which point the posi- 
tion and momentum 4-vectors are saved. All electrons 
are tracked until their energy decreases to less than 75 
GeV, or until they reach the surface of z=0. Particular 
care is given to the computation of time delays. The 
time delay of individual particles is computed relative 
to the arrival time of a putative photon propagating di- 
rectly from the source to the current position of the par- 
ticle. This procedure is adopted to maintain precision in 
numerical simulations down to minute timescales. The 
computation accuracy of time delays was verified using 
an arbritrary precision numerical integrator developed 
at Lawrence Berkeley National Laboratory The full 
solution for the equations of motion of electrons in a cos- 

1 http: //crd- legacy . lbl .gov/- dhbailey/mpdist/ 
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Figure 5. Skymaps for source at z=0.13, a = 1.5, E c = 30 TeV, and B IG mf=10 15 , ohs = 0° (upper left), 2° (upper right), 5° (lower 
left), 10° (lower right) 



mologically expanding universe and under the influence 
of constant magnetic field are used to track the position 
and momentum 4- vectors between IC interactions. 

Previous models ranging from analytic to Monte 
Carlo codes have been reported in the literature to 
genera t e constraints on the I G MF dNeronoy fc V ovk 
(12Q1QD; Ij avecchi o et all (12010ft: iDermer et all (12011ft: 
Esse v et all (12010ft : Dolag efal] (|2011ft : iTavlor et all 
(|2011ft : iHuan et all (|2011h ). Each one has utilized vari- 
ous degrees of simplification through the use of solutions 
of one or two dimensional kinetic equations or simpli- 
fied analytical approximations of QED interactions or 
geometrical effects of cascade development. For exam- 
ple, nearly all of these approaches neglected the Klein- 
Nishina regime of IC scattering and ignored the EBL as a 
target photon field for this process. The VHE secondary 
photons produced by the highest energy electrons are ca- 
pable of generating second and higher orders of the EM 
cascades, and this effect is neglected in most of the prior 



works. 

Figure [6] illustrates the importance of higher order cas- 
cading to correctly describe secondary photons with en- 
ergies > 200 GeV. This figure shows the results of the 
simulations of secondary emission produced for a source 
at z=0.3, with B IGM f 10" 16 G, V 0, and T = 10. 
The two panels display different intrinsic spectra for the 
source; panel a) is obtained with a = 1.5 and e c = 30 
TeV, and panel b) corresponds to a = 1.5 and e c = 5 
TeV. Both panels show the direct differential flux energy 
density (dFED) of the source together with the time- 
integrated secondary dFED, within 0.5° from the posi- 
tion of the source. The two lines of the secondary dFED 
shown on the figure are obtained with second order cas- 
cading included or not included in the simulation. It is 
evident that the secondary dFED > 200 GeV is strongly 
affected by this assumption and if this simplification is 
made, the Monte Carlo simulations may over-predict the 
total dFED in the VHE energy range. 
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Figure 6. Comparison of spectra produced by both neglecting and including multiple generations in the cascade for a source at z=0.3, 
with Bjgmf = le-16 G, # v = 0, and F = 10. (integrated over a 0.5 deg aperture). The black (solid) line is the prompt, direct spectrum, 
the red (long dashed) line shows the secondary spectrum produced when multiple generations of the cascade are included, and the green 
(small dashed) line shows the secondary spectrum when only a single generation of cascading is allowed in the simulations. Two different 
source spectra were used a) (left) spectral index a = 1.5, cutoff energy, e c = 30 TeV and b) (right) a = 1.5, e c = 5 TeV. 



3. GAMMA RAY DATA, INSTRUMENTATION, AND 
STRATEGY FOR DATA ANALYSIS 

The gamma-ray data used in this study are compiled 
from the data sets of two types of instruments. In the 
HE regime, the data were obtained by the Fermi-L AT in- 
strument and the data were processed utilizing publically 
available tools, version v9r23pl, with the update from 
November 6, 2011. For the VHE regime, data previously 
reported by the Imaging Atmospheric Cherenkov Tele- 
scopes (IACT) instruments VERITAS and HESS were 
used. In both the HE and VHE energy regimes, the 
software used for analysis interprets the gamma-ray flux 
as originating from a point source; no angular extension 
is assumed. Since the secondary photons of intergalac- 
tic cascades inherently represent an extended source, 
the methods of comparison of simulations and data are 
instrument-specific. 

3.1. VHE Data Considerations 

In most of the following discussion, the validity of the 
Bigmf = hypothesis (hereafter called HO) is exam- 
ined. Under this assumption, the source of the gamma- 
ray photons has an angular extent due to the QED pair 
production and IC scattering processes. The angular 
size of this extension is much smaller than the gamma- 
ray point spread function (PSF) of both IACTs and the 
Fermi-L AT instrument, and therefore, the gamma-ray 
sources are point-source like. In the VHE regime, the 
point source assumption becomes invalid for Bigmf ^ 
1CT 15 G because the PSF of IACTs and the angular ex- 
tension of the source become comparable. This require- 
ment, however, is dependent on the higher energy end of 
the spectral energy density (> 10 TeV) of the source and 
its history of activity. 

To process simulated data in the VHE energy regime, 
the instrumental PSF of 

is used as suggested in lAharonian et al.l (|2006d ). The 
typical values of parameters of the PSF are 9\ = 



0.06505°, 6> 2 = 0.1697°, a = 0.505. These parameters are 
assumed to be energy independent in analyses of IACT 
data reported and therefore, we make a similar assump- 
tion for processing simulated data. For each simulated 
photon, Equation[3]is used to find the probability for this 
photon to be reconstructed within the 68 % containment 
radius from the position of the putative point source. 
The effective point source flux from an AGN is estimated 
as the total simulated flux within the 68 % containment 
radius divided by 0.68. This method of flux evaluation 
for each energy bin accurately models the IACT data as 
reported in the literature. 

The VHE data set used in this paper is summarized in 
Table [TJ The data sets of the first three sources (RGB 
J0710+591, 1ES 121 8+304, and 1ES 022 9+200) are iden- 
tical to those used in iTavlor et al.l (|2011[ ) , except that an 
additional data set for 1ES 1218+304 obtained by the 
VERITAS Collaboration just before the launch of the 
Fermi satellite (on August 4, 2008 or M JD 54682) is 
conside red in this study. As reported in lAcciari et al.l 
(|2010aj ). the activity of the source is nearly identical 
during these non-overlapping periods, except for an el- 
evated flux of the source peaking at the level of ~20 
% Crab over a few nights of observations. The data 
set for 1ES 0229+200 was also obtained prior to the 
launch of the Fermi satellite. Based on th e report from 
iPerkins fc VERITAS Collaboration! (|2Q10f ). the activity 
of the source as measured by VERITAS during the sec- 
ond year of the Fermi mission, appeared to resemble 
the reported SED by the HE SS collaboration prior t o 
the launch of Fermi satellite (jAharonian et al.l l20Q7ch . 
VERITAS has continued monitoring this source since the 
Fermi launch and tentatively detected flux variations on 
a yearly time scale (private communication, J.S. Perkins 
and VERITAS collaboration). Finally the data sets for 
four other extreme TeV blazars (1ES 0347-121, 1ES 1101- 
232, H 2356-309, and RGB J0152+017) were taken from 
the discovery publications by the HESS collaboration, 
and all of these sources were observed prior to the start 
of the Fermi mission. 

3.2. HE Data Considerations 
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Table 1 

Data of RGB J0710+591 are taken by VERITAS TAcciari et all (I2010bl) Source 1ES 1218+304 wa s observed by 
VERITAS and the results were published in two separate submissions Acciari et al. ( 2010a, 20 Q3). This sourc e 
was also detected during a six day period around MJD 53750 by the MAGIC collaboration I Albert et aLl (2006), 
at a level similar to that of the reported flux by VERITAS about one year later. Observations of 1ES 0229+200, 
1ES 0347-121, 1ES 1101-232, H 2356-309, and RGB J0152 +017 were taken prior to the launch of the Fermi 
satellite, and are reported by the HESS collaboration Aharonian et al. (2007c. b. a, 2006b, 2008). 



IACT data summary 


Source 


z 


IACT 


Flux [10- 12 cm-^ 


1 s- 1 ] 


index 


MJD (approx) 


Hrs 


RGB J0710+591 


0.125 


VERITAS 


(> 300 GeV) 3.9 


± 0.8 


2.7 


± 


0.3 


54882-54892 


22.1 


1ES 1218+304 


0.182 


VERITAS 


(> 200 GeV) 18.4 


± 0.9 


3.1 


± 


0.3 


54829-54944 


27.2 


1ES 1218+304 


0.182 


VERITAS 


(> 200 GeV) 12.2 


± 2.6 


3.1 


± 


0.1 


54070-54220 


17.4 


1ES 0229+200 


0.14 


HESS 


(> 580 GeV) 0.94 


± 0.24 


2.5 


± 


0.2 


53614-53649, 53967-54088 


41.8 


1ES 0347-121 


0.188 


HESS 


(> 250 GeV) 3.91 


± 1.1 


3.1 


± 


0.3 


53948-54100 


25 


1ES 1101-232 


0.186 


HESS 


(> 200 GeV) 4.5 


± 1.2 


2.9 


± 


0.2 


53111-53445 


43 


H 2356-309 


0.165 


HESS 


(> 200 GeV) 4.1 


± 0.5 


3.1 


± 


0.3 


53157-53370 


40 


RBG J0152+017 


0.08 


HESS 


(> 300 GeV) 2.7 


± 1.0 


2.9 


± 


0.5 


54403-54418 


15 



For an AGN with significant power emitted at E > 
few TeV the angular extent of the source due to cas- 
cade emission may become comparable to the PSF of 
the Fermi-L AT at fields as low as Bigmf ^ 10 -17 G. 
For Bigmf magnitudes less than this, an AGN is effec- 
tively a point source for the LAT. To evaluate the effec- 
tive point source flux of an AGN for the Fermi-L AT, a 
procedure similar to that of the IACT case was adopted. 
For every simulated photon of a given energy, the energy- 
dependent PSF of the Fermi-L AT was used to determine 
the probability of reconstruction of this photon within 
the 68 % containment radius from the position of the 
putative point source. The flux evaluated within the 68 
% containment was again rescaled by 0.68 -1 to estimate 
the effective point source flux in each energy bin. The 
PSF used for this conversion was determined from a 2 
year time-averaged sample of AGIS0. 

To compare simulated differential fluxes from an effec- 
tive point source to the Fermi-LAT data, the standard 
analysis tools were applied but with a notable important 
distinction from previous studies, which are cited in sec- 
tion [2]4l Since in the HE regime, the flux of gamma-ray 
photons can be dominated by either prompt or secondary 
emission, we first derive from simulations, the spectral in- 
dex in each energy bin. This index is then used as a fixed 
parameter for the maximal likelihood evaluation^] of the 
flux in each energy bin in the Fermi data, within the 10° 
region of interest (ROI) which also includes all nearby 
sources from the 2 year point source catalog and diffuse 
backgrounds. The HE point source fluxes or upper limits 
are derived using this procedure. 

3.3. Data and Model Comparison 

To compare the predictions of the Monte Carlo simula- 
tions to the data, we combined a x 2 -hke parameter, that 
includes both the HE and VHE regimes. The simulated 
effective point-source flux is used as the model expecta- 
tion value, and we assume that the statistical error based 
on a vast amount of simulations is negligible compared 
to the observational error. The point-source fluxes de- 
rived with the use of the Fermi tools as explained above 

2 http : //f ermi . gsf c . nas a. gov/ ssc/dat a/ analysis/document at ion/ 
|Pass7_usage . html | 

a http: //f ermi .gsf c .nasa.gov/ssc/data/analysis/ 



or obtained from IACT publications are used as data. 
It is assumed that the observational error obtained or 
reported is the primary source of discrepancy between 
the data and the model. A x 2 -like parameter is used to 
estimate goodness-of-fit and convert it to a confidence 
level (or the probability of exclusion) assuming that the 
errors in both the HE and VHE regimes are dominated 
by statistics. Throughout the paper, we test the null hy- 
pothesis, HO, that Bigmf — is incompatible with the 
data, and the confidence level is the probability that the 
result cannot be obtained with this assumption. 

The model spectral energy density is derived based 
on the full Monte Carlo simulations computed for mo- 
noenergetic primary photons with 8 bins per decade, of 
equal width in logarithmic space. The simulated spectral 
energy density data are equally binned with 8 bins per 
decade and each bin centered on the energy of the pri- 
mary photon monoenergetic line. The VHE data are re- 
ported in different publications at different energies and 
with different binning. We use simulated data to inter- 
polate the flux value to the reported positions of the bins 
and their widths. In the 3 decades of the HE regime, 4-6 
energy bins are generated depending on the given source 
luminosity and statistics. The simulated data are then 
used to interpolate the flux value to these energy points. 

To compare HE and VHE data of each source with the 
simulations, effective point source fluxes are generated 
for a set of gamma-ray source models with fixed values 
of a, e c , T, and V . For each model, a is chosen in the 
range (1.0, 2.5) with a step of 0.1, and e c is chosen in the 
range (600 GeV, 60 TeV) with 8 bins per decade, equally 
spaced in logarithmic energy. Each model is character- 
ized with default values of V = 0, and T = 10. These 
two parameters are varied only when systematic effects 
of an individual source are explored. The flux normal- 
ization factor, Fo, is determined by minimization of the 
VHE part of x 2 by solving dxy HE /dF$ = 0. The spec- 
tral break energy e#, and the spectral index 7, relevant 
to the HE part of the spectrum are chosen by minimizing 
the overall x 2 value, by allowing cb to vary between >10 
GeV to the lower edge of the lowest energy VHE data bin 
and 7 between -5 and 5. The x 2 value for each model 
was converted into a confidence level, using x 2 statistics 
with N - 3 degrees of freedom, where N is the number of 
data points in both the HE and VHE regime, and there 
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are three parameters fit for each model (Fo, e#, and 7). 

As an illustration of a typical result of the data and 
model comparison, Figure [7^ shows the \ 2 confidence 
level of the dFED. For this set of models with a range of 
a and e c , the most favored model has a value of a = 1.8 
and e c = 3.16 TeV and is incompatible with the data at 
the level of about 75 %. The computations of x 2 for this 
fit assumed HO, and Xigmf — 1 Mpc. Figure [Tb shows 
the data points for the dFED and the best fit simula- 
tion result with these a, e c parameters. For this source's 
dFED, below 10 GeV, the spectrum is dominated by cas- 
cade emission, while above 10 GeV, it is dominated by 
prompt radiation. 

4. IGMF CONSTRAINTS AND EFFECTS OF SYSTEMATIC 
UNCERTAINTIES 

The primary goal of this paper is to pro- 
vide precise numerical verification o f the constraints 

on th e IGMF as rep o rted in Nerono v fc Voykl 

(f20T0h: iTavecchio et al.l (l20Toh: Ipolag et al.l (|20ll ): 
iDermer et al.l ([201 lh : iTavlor et al.l ([201 If ), since most of 
these results were obtained with various simplifications in 
the analysis, and some were derived with semi- analytical 
approach es to qualitat i vely v erify more detailed compu- 
tations ([Dermer et al.l ([20 111 )). Particular focus is given 
to establishing the robustness of magnetic field limits 
when various systematics are taken into consideration. 

Three sources (RGB J0710+591, 1ES 1218+304, and 
1ES 0229+200) are used in this study which have been 
reported as having provided c onstraints on the I GMF, 
in the comprehensive study of ITavlor et al.l (|2Qllh (and 
references within). Hereafter, we refer to this paper as 
TVN11. An additional 4 hard-spectrum TeV blazars ob- 
served prior to the beginning of the Fermi mission are 
considered-lES 0347-121, 1ES 1101-232, H 2356-309, and 
RGB J0152+017. The Fermi-LAT data for these sources 
are re-analyzed and are collected from the mission start 
time Aug 4, 2008 to February 14, 2012, and the updated 
P7SOURCE IRFs are used along with the Pass 7 data. 
The models for extragalactic and diffuse backgrounds 
were used together with the standard gamma-ray selec- 
tion constraint of zenith angle < 100° which eliminates 
earth limb gamma-rays. 

Perhaps the main source of uncertainty in placing con- 
straints on the IGMF stems from the unknown duty cy- 
cle of TeV blazars and particularly, the history of the 
hi ghest energy TeV em ission, as has been pointed out 
in IDermer et al.l ([201 lh . The sampling of the VHE ac- 
tivity of these sources reported by IACTs is limited to a 
few tens of hours dispersed over a period of a few weeks 
to a few years. In the regime of very low IGMF (Bigmf 
< 10 -20 G), most of the secondary radiation from inter- 
galactic cascades with energy > 100 MeV, which orig- 
inates from the primary VHE flux sampled by IACTs, 
would have reached the earth and would be detected by 
the Fermi-LAT (see Fig. |2j) within a few hours. This 
assumes that the HE flux from a given source sampled 
by the Fermi-LAT over the period of the mission (about 
4 years) could be viewed as "contemporaneous" to IACT 
measurements, for the purposes of verification of HO. 
This explicitly assumes, however, that the duty cycle of 
a given source in the VHE regime = 1 over this same 
period. 



4.1. Analysis of RGB J0710+591 Data 

The VHE observations of RBG J0710+591 are sum- 
marized in Table [TJ They include 5 energ y data points 
reporte d by the VERITAS collaboration in lAcciari et al.l 
(|2010b[ ), observed during the time period December 2008 
- March 2009 for a total of 22 hours. For the HE regime, 
6 energy bins were used spanning from 200 MeV to 200 
GeV. A detection in each energy bin was assumed if the 
Test Statistic (TS) > 9 with the fixed value of spectral 
index in the bin, predicted from simulations. Otherwise, 
an upper limit was computed. It is important to note 
that with this strategy and extended data set as com- 
pared to previously used in TVN11, the flux in the low- 
est energy bin now constitutes a detection rather than 
an upper limit. This data point had been critical for 
rejecting HO. 

Figure^ shows the confidence level, for HO, for a set of 
models c harac terized by the range of e c and a described 
in section [3731 It identifies the best fit model with values 
of a = 1.8 and e c = 3.16 TeV, which is incompatible 
with HO at the < 80 % confidence level. The fit of the 
simulated dFED and observations is shown in Figure [7]o. 
It appears that the conclusion of TVN11 that HO is ruled 
out at the 98.8 % level is invalidated due to two factors. 
First, the Fermi-LAT dataset underwent revision from 
the old pass 6 version (P6) to the current pass 7 (P7) 
and this allowed a detection to be made in the lowest 
energy bin, which is higher than the previously computed 
upper limit. The second factor may be due to the fitting 
algorithm which is different from that adopted in TVN11, 
in which the flux or upper limit determination in a given 
energy bin was fixed to the best fit index over the entire 
energy range. 

Furthermore, the geometrical orientation of the jet 
with respect to the observer and the jet boost factor rep- 
resents another source of uncertainty, and tuning these 
parameters can further improve the goodness of the x 2 
fit. For example, TVN11 assumes a viewing angle of 2° 
and an effective jet opening angle of 6°, corresponding 
to a boost factor of ~ 10. As shown in Figure HJ how- 
ever, lower boost factors or smaller viewing angles lead 
to lower total power of the jet at the highest energies, 
and therefore lead to reduced secondary flux. 

4.2. Analysis of 1ES 1218+304 Data 

The VHE observations of 1ES 1218+304 are summa- 
rized in Table [TJ which includes 2 data sets. The first set, 
obtained during December 2008 - May 2009, is based on 
27 hours of data and has 9 energ y data points r eporte d 
by the VERITAS collaboration in lAcciari etHI (j2010at ). 
This data set was used in TVN11, and for consistency 
it is also used in this study. The Fermi-LAT data for 
this source were produced in the same way as for RGB 
J0710+591, with 6 energy bins spanning from 200 MeV 
to 200 GeV. This source has excellent statistics in each 
Fermi-LAT energy bin, with TS > 25. It is important 
to note that the extended exposure and updated pass 7 
data set used in this work shows no statistically signifi- 
cant difference compared to that of TVN11. 

Figure [8^ shows the confidence level on the a-e c plane, 
obtained for HO. It suggests a best fit model with values 
of a = 1.8, e c = 3.16 TeV, a spectral break energy of €b 
= 10 GeV, and an index below the break energy of 7 = 
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Figure 7. Source RGB J0710+591 analyzed under the assumption of HO. (a) (left) Confidence level for a given a, and ec, where the 
remaining three parameters of the gamma-ray spectral model (see Eq. U-7, es, and Fo) are chosen so as to minimize % 2 . (b) (right) 
Simulated dFED of the best fit model, of a = 1.8, e c = 3.16 TeV, showing both the prompt and secondary cascade contributions to the 
total dFED, along with the HE and VHE data. 

1.0, which is incompatible with HO at the less than 70 
% confidence level. The fit of the simulated dFED and 
observations is shown in Figure [SJd. It is evident that the 
conclusion of TVN11, that HO is ruled out with more 
than 99.99 % probability, is purely due to the assumption 
that a single source spectral index holds for over five 
orders of magnitude in energy. Allowing a spectral break 
energy and an intrinsic spectral index below the break 
energy to vary as detailed in section [231 makes it possible 
to interpret the 1ES 1218+304 data set as consistent with 
HO. 

The amount of secondary radiation strongly depends 
on the power output of the TeV blazar at the highest 
energies. For this source, more so than the others, IACT 
observations demonstrate strong variability in the VHE 
spectrum. The VERITAS collaboration reports that the 
1ES 1218+304 data were sampled sparsely over a period 
of 115 days in late 2008 - 2009, and while the majority 
of the data are consistent with a steady baseline flux, 
the data set also includes a statistically significant flare 
which peaked at ~ 20 % Crab, and lasted a few nights. 
The flux at the peak of the flare was 3-4 times higher 
than the baseline flux and it significantly increases the 
average flux value observed over the entire period. Fur- 
thermore, evidence for variability of this source can be 
inferred from the VERITAS publication covering its 2 
year activity which o ccurred prior to the Fermi mission 
(jAcciari et all ({2009)). The flux observed at that time 
constitutes about 60-70 % of that reported in the second 
data set (see Table [T]). Overall, the IACT data to date 
suggest that the average observed VHE flux of this souce 
could be lower than used in TVN11, yet the assumption 
of the higher average VHE flux is still compatible with 
HO. 

4.3. Analysis of 1ES 0229+200 Data 

The parameters of the data set of 1ES 0229+200 are 
given in Table Q] and include two sets of observations by 
the HESS collaboration during the period from Septem- 
ber 1, 2005 to December 19, 2006, acc umulating 41.8 
hours exposure (jAharonian et al.l (j2007cf )). This data set 
provides the time-averaged dFED for 8 bins over the en- 



ergy range spanning from 500 GeV to 16 TeV. This same 
data set was used in the previous study of TVN11. The 
Fermi-LAT dFED for 1ES 0229+200 utilizes four evenly 
spaced bins in log space in the range from 420 MeV to 300 
GeV. Only the first energy bin in this data set provides 
a strong detection (TS > 25), for all simulated models in 
which the secondary flux dominates the total flux. The 
TS for all other energy bins is typically found at > 9 
for the majority of simulated models but in some cases, 
only the upper limit can be established (TS < 9). This 
indicates that the source detection in the HE regime is 
weak. Perhaps more so than for any other source, the 
dFED of 1ES 0229+200 does not resemble a power law 
in the HE energy regime, making the y 2 fits relatively 
poor. 

To evaluate the goodness of fit of the data to the Monte 
Carlo simulations, we assume that the data accumulated 
by the HESS collaboration over 2005-2006 is representa- 
tive of the source activity during the first 3.5 years of the 
Fermi-LAT data used in this work. The x 2 fit obtained 
under this assumption and for HO is shown in Figure [9^. 
We confirm the result of TVN11 that this source does 
not have a viable source model that explains the com- 
bined HE- VHE data set and agree that HO is ruled out 
at the 99.5 % confidence level. The best fit (a = 1.3, 
e c = 1 TeV) model requires a dramatic spectral break 
just below 100 GeV and the dFED of 1ES 0229+200 be- 
low this energy is completely dominated by secondary 
flux as shown in Figure [9]o. 

The effects of two systematic uncertainties, viewing an- 
gle V and Doppler factor T, were considered. The de- 
fault assumptions in the analysis are 6 V = and T = 10. 
Increasing the viewing angle, e.g. 6 V = 2° as used in 
TVN11, would further overproduce radiation in the HE 
regime (see Figure H^i). Increasing the Doppler factor, T, 
combined with the V = assumption would imply that 
the overall luminosity of the jet in the VHE band should 
be lower to fit the VHE observations since the jet is colli- 
mated into a smaller angle. Rescaling of the prompt ra- 
diation to fit the VHE data will however equally rescale 
secondary emission in the HE regime if the angular distri- 
bution of the prompt photons is significantly wider than 
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Figure 8. Source 1ES 1218+304 analyzed under the assumption of HO. (a) (left) confidence level for a given a, and e c (b) (right) 
simulated dFED of the best fit model, of a = 1.8, e c = 3.16 TeV, showing both the prompt and secondary cascade contributions to the 
total dFED, along with the HE and VHE data. 
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Figure 9. Source 1ES 0229+200 analyzed under the assumption of HO. (a) (left) confidence level for a given a, and e c . (b) (right) 
Simulated dFED of the best fit model (C.L. = 0.995) of a = 1.3, e c = 1 TeV, snowing both the prompt and secondary cascade contributions 
to the total dFED, along with the HE and VHE data. 



the characteristic scattering angles acquired in the QED 
processes. Therefore, for reasonable values of T < 100, 
the change in the secondary radiation above 100 MeV 
for the best fit model was found to be negligible. The 
effect becomes of order 10 % in the lower part of the HE 
spectral range only for T ~ 10 4 - 10 5 , and it cannot be 
used to reconcile the observational data of 1ES 0229+200 
with HO. 

Since the energy density of the EBL directly affects 
the propagation length of VHE photons, the uncertainty 
in the EBL model represents perhaps the most impor- 
tant source of systematic error. To investigate this, 
two additional EBL models were generated, EBL model 
2 & 3, shown in Figure [TJ EBL model 2 is charac- 
terized by a considerably lower energy density in the 
far infrared peak of the dust emission. This model 
is motivated by the recently resolved lower limit on 
the EBL which is based on t he galaxy counts in th e 
data obtained wit h the Spitzer jBethermin et al. (2010); 
iDole et all (120061)). Her sc hel dBerta et al.l (120101)). and 
AKARI (jMatsuura et all (|2011|)) satellites. This model 



corresponds to the lowest possible far infrared EBL en- 
ergy density allowed within 2a. It was found that even 
with such an extreme assumption about the far infrared 
EBL, the decrease of the secondary radiation in the 
model of 1ES 0229+200 was negligible. This conclu- 
sion is due to the fact that the source models providing 
the best x 2 fits have high energy cutoffs, e c , below ~ 5 
TeV, and are thus insensitive to EBL photon wavelengths 
> 25/im (kinematic threshold of pair production). 

The EBL model 3 is based on the resolved EBL 
energy density of the starlight peak, which is de- 
rived from galaxy c ounts utilizing data from the HST 
(jMadau fc Pozzettil (2000)) in the visible ( from 0.36 



2.2 £m i), Spitzer in the near-IR (3.6 - 8 /im) (I Fazio et al 
dlOOl) and ISO in the mid-I R (15-24 /im) (lElbaz et al 



( 20020: IPapovich et al.l ([20041 )). As compared to the de- 
fault model 1, this EBL model has an energy density in 
the visible reduced by about 25% and an energy density 
in the near-IR reduced by about 50%. The model is com- 
patible with the galaxy counts results to within la but it 
effectively does not allow any additional contribution to 
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Figure 10. Source 1ES 0229+200 analyzed, assuming EBL Model 3 and HO. (a) (left) Confidence level for a given a, and e c . (b) 
(right) Simulated dFED of the best fit model (C.L. = 0.78) of a = 1.6, e c = 3.16 TeV, showing both the prompt and secondary cascade 
contributions to the total dFED, along with the HE and VHE data. 

fact, this feature was used to derive upper limits on the 
EBL energy density in the mid-IR, which are taken to 
be the values of model 1, by assuming that the spec- 
tral index of t he source 1ES 11 1-232 cannot be harder 
than 1.5 (e.g. lAharonian et al.l (|2006a[ )). The 25 - 50% 
lower mid-IR density of model 2 significantly softens the 
intrinsic spectrum of 1ES 0229+200 reducing the total 
energy available for the development of the intergalactic 
cascade, and therefore the flux of the secondary photons. 
Thus, there are a range of EBL models with mid-IR en- 
ergy density bounded by the lower limits on the EBL 
to some SED slightly below that of model 1 which are 
compatible with the EBL lower limits and HO. 

In a recent study of the dual constraints on the EBL 
and IGMF it was found that an EBL model similar to 
model 3, is still incompatible with HO and w ould require 
a lowe r bound of Bjgmf = 6 x 10 -18 G (|Vovk et al.l 
(2012)). All source models analyzed in that paper had 
a single cutoff energy ec = 5 TeV and varying spectral 
index a in the range of - 1.5 with a single power law 
dFED over the entire VHE and HE energy range. A 
similar assumption of a steady flux from 1ES 0229+200 
over the lifetime of the Fermi-L AT was made. Figure 
ITT1 illustrates the confidence level for a wider range of 
source models analyzed in this work together wi th the 
range of models considered bv IVovk et al.l (j2012[ ). The 
figure confirms the incompatibility of HO with the data 
given assumptions about the source model used in that 
paper. However, in an extended parameter space of the 
1ES 0229+200 models, HO can be reconciled with obser- 
vations of this source. 

Another avenue to make the HE- VHE observational 
data of 1ES 0229+200 compatible with HO is to question 
whether or not the VHE flux of the source reported is rep- 
resentative of the average flux during the Fermi mission. 
The HESS measurements were accumulated in 2005 (6.8 
hrs) - 2006 (35 hrs) and are not strictly contemporaneous 
with the Fermi HE spectrum. The significance of the de- 
tection in 2005 reported is 2.7<r, while in 2006 it is 6.1<j 
with average photon fluxes above 580 GeV of 6.8 x 10 -13 



Figure 11. Confidence level for the simulated source models 
for 1ES 0229+200 obtained under the assumption of EBL Model 
3. The black line t erminated with cr osses represents the range of 
models analyzed in Vovk et al. (20l3)- 

the EBL from unresolved or unknown sources. The re- 
sults of the simulations of intergalactic cascading for this 
model with HO are shown in Figure [TOh . It was found 
that a number of 1ES 0229+200 source models are com- 
patible with the combined VHE and HE data set, and 
one of the best fit examples characterized by a = 1.6, 
e c = 3.16 TeV is shown in Figure [TUb. 

The strong sensitivity of the secondary photon flux to 
the EBL energy density in the near-IR combined with 
the conspicuous lack of a non-trivial EBL absorption fea- 
ture in the VHE energy band (~ 200 GeV - 5 TeV) is 
due to the peculiar behavior of the EBL in this wave- 
length range. For XI \ oc A -1 , the optical depth is inde- 
pendent of the energy of a VHE photon (gray opacity). 
A small deviation from this proportionality results in a 
logarithmically slow dependence of the optical depth on 
the energy of the VHE photon, producing a power law, 
rather than expo nential-li k e ch ange in a blazar spectrum 
as discussed in Vassiliev (20(X)j). The behavior of the 
SED in the mid-IR exactly satisfies this condition and 
explains the "invisibility" of EBL absorption effects in a 
blazar dFED. The effect however is strong and reflected 
in the change of the spectral index of this blazar. In 



cm -2 s _1 , and 10 xlO -13 cm -2 s _1 , respectively. Due to 
the low flux of this source (1-2% of the Crab nebula flux), 
and the small data set in the original 1ES 0229+200 dis- 
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Figure 12. Spectral model of 1ES 0229+200 of a = 1.3, e c = 1.78 
TeV (C.L. < 0.8), for the TeV end of the spectrum scaled down 
by 10~ 1//2 of their flux values to that reported in Aharonian et all 
(20073). 

covery paper, statistically significant flux variability as 
observed in 2005 and 2006 was not detected. Although 
these observations are compatible with a constant flux, 
the variability hypothesis cannot be ruled out, based on 
the statistical and systematic errors reported. Further- 
more, obser vations of this source in 2009, as reported b y 
VERITAS (jPerkins fc VERITAS Collaboration! (pQlOh ). 
were compatible with the average flux value of the HESS 
data set. However, the average flux obtained was dom- 
inated by a period of significantly higher "flaring" ac- 
tivity during a single dark run. In general, VHE obser- 
vations above a few TeV (relevant for secondary photon 
production) require considerable integration time and so 
far, they are too sparse to claim that the HESS value of 
the flux is representative of the average flux during the 
Fermi mission. Further communication with the VER- 
ITAS Collaboration suggests that the flux level of 1ES 
0229+200 has been steadily declining from 2009 - 2012 
(private communication). To investigate the effect of a 
reduced duty cycle for 1ES 0229+200, the spectrum of 
this source was modified at the highest 5 energy points 
to half an order of magnitude of their reported values. It 
was found that the VHE - HE data set combined in this 
way does not rule out HO at more than 95% confidence 
level. One of these compatible models, with a = 1.3 e c 
= 1.78 is shown in Figure fT2l Therefore, the conclusion 
that the HO is ruled out with high significance heavily 
rests on the assumption that the HESS measurements are 
representative of the average flux for E > 2 TeV and a ~ 
10 _1//2 change of the highest energy part of the spectrum 
invalidates this conclusion. 

4.4. Analysis of 1ES 0347-121, 1ES 1101-232, H 
2356-309, and RGB J0152+017 

Prev i ously published IGMF studies dNeronov fc Vovkl 
(|2010[ ): lEssev et al.l (|2010f ): iTavecchio et al.l (|20f!f )T have 
used additional extreme TeV blazars to derive con- 
straints on the IGMF, four of which are considered in 
this section, 1ES 0347-121, 1ES 1101-232, H 2356-309, 
and RGB J0152+017. All of these sources were observed 
by the HESS collaboration prior to the beginning of the 
Fermi mission. Therefore, to investigate HO, it is neces- 



sary to assume that the VHE activity of these sources as 
characterized by the HESS collaboration typically during 
the period of 2004 - 2007 is representative of the VHE 
activity over the duration of the Fermi mission. The pa- 
rameters of these data sets are summarized in Table [TJ 
The Fermi-LAT time-averaged dFED for all sources was 
derived from the beginning of the mission until February 
14, 2012, and depending on the strength of the source, 
was computed for either 4 or 6 evenly spaced logarithmic 
bins over the energy range of 200 MeV - 200 GeV. 

The VHE data set of 1ES 0347-121 consists of a set of 
observations over the period of August - December 2006, 
during which time 25.4 hours exposure was acculumlated. 
A time-averaged dFED for 7 bins over the energy range 
from 25 GeV - 3.67 TeV was derived (I Aharonian et al.l 
(|2007bf )V The flux found in the first and fourth (last) 
bins is weak (TS < 9) for all spectral indices tested, al- 
lowing only an upper limit to be established. The flux 
in the second and third bins is typically found with 9 < 
TS < 25 for the simulated models where the secondary 
flux dominates the total flux. Figure [T3a shows the con- 
fidence level of simulated models obtained for the HO. 
The best fit models in the a - ec plane are found at ec 
values near 1 TeV, and the dFED for one of these models 
is illustrated in Figure [T3b. The relatively low confidence 
level of the 1ES 0347-121 simulated models is partially 
due to the poor fit of the highest energy bins of the VHE 
regime where the reported dFED tentatively exhibits a 
feature of increasing energy density. This trend in the 
dFED is not accounted for in the set of simulated mod- 
els investigated. A similar spectral feature appears to 
be even more pronounced in the VHE data set of 1ES 
1101-232, which perhaps may signal unmodeled physics 
process(es) or a systematic error in the data analyses. 

The VHE data set of 1ES 1101-232 consists of 3 pe- 
riods of observations spanning from April 2004 - March 
2005, for a total of 43 hours. The time- averaged dFED 
is reported for 10 bins over the energy range 225 GeV - 
4 TeV. The Fermi-LAT dFED for 1ES 1101-232 is espe- 
cially weak, and in fact, the first bin allows a determi- 
nation of only an upper limit (TS < 9) for all simulated 
models tested. Figures [Till and b illustrate the compat- 
ibility of the 1ES 1101-232 VHE and HE data sets with 
HO, despite the fact that the previously described fea- 
ture in the highest energy bins of the VHE regime is not 
well fitted by the models. 

Observations of H 2356-309 were obtained over the pe- 
riod of June - September 2004, for a total exposure of 40 
hours. The time-averaged dFED was provided for eight 
bins over the energy ran ge from 200 GeV - 1.23 TeV 
(j Aharonian et al.l (j200 6b)). The flux in the first energy 
bin is weakly detected for most simulated models (9 < 
TS < 25), but the second and third bins exhibit a strong 
detection (TS > 25). An upper limit is derived for the 
fourth bin in the HE dFED due to a weak signal present 
(TS < 9). As illustrated in Figure [T5| this source has a 
very large set of models compatible with the HO. Most 
of these models, however, suggest that the flux in the 
two lowest energy bins of the HE regime is dominated by 
secondary radiation. 

RGB J0152+017 was observed over the period of Oc- 
tober 30 - November 14 2007 for a total exposure of 
14.7 hours. A time averaged dFED was reported for 
six bins over the energy range of 240 GeV - 3.6 TeV 
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Figure 13. Source 1ES 0347-121 analyzed under the assumption of HO. (a) (left) Confidence level for a given a, and e c . (b) (right) 
Simulated dFED of the best fit model of a = 1.3, e c = 0.75 TeV, showing both the prompt and secondary cascade contributions to the 
total dFED, along with the HE and VHE data. 




Figure 14. Source 1ES 1101-232 analyzed under the assumption of HO. (a) (left) Confidence level for a given a, and e c . (b) (right) 
Simulated dFED of the best fit model of a = 1.2, e c = 0.75 TeV, showing both the prompt and secondary cascade contributions to the 
total dFED, along with the HE and VHE data. 




Figure 15. Source H 2356-309 analyzed under the assumption of HO. (a) (left) Confidence level for a given a, and e c . (b) (right) 
Simulated dFED of the best fit model of a = 1.8, e c = 4.22 TeV, showing both the prompt and secondary cascade contributions to the 
total dFED, along with the HE and VHE data. 
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(jAharonian et al.l (|2006bl )). The flux in all but the first 
energy bin of the Fermi data is strongly detected (TS > 
25) for the majority of simulated models. Figure [T6l a and 
b illustrate that this source perhaps provides the weakest 
constraints on the IGMF with nearly the entire parame- 
ter space of simulated models compatible with the VHE 
and HE spectral data. In only some of these models, the 
HE part of the spectrum is dominated by the primary 
emission. 

5. DISCUSSION 

In this paper, we have investigated the HE - VHE en- 
ergy spectrum of seven extreme TeV blazars for which 
radiation in the HE band may be dominated by the 
secondary photons produced through cascading. These 
sources are characterized by their observed hard spectra 
in the VHE band accompanied by redshifts of order ~ 
0.1 suggesting a very large energy output into pair pro- 
duction and subsequent cascading. For these sources, 
observations in the HE band can therefore limit the flux 
of secondary photons and establish a lower bound on the 
IGMF . This strategy ha s been u tilized in several publica - 
tions dNeronov fc Vovlj (|2Q1QD ; [Tavecchio et al] (|2Q1Q|); 
iDermer et al.l (120111 ): lEssev et al.l (l2010f): iDolag et al.l 
(|201ir >: lTavlor et al.l (j2Qllh : iHuan et al.l (|2QllD ) to sug- 
gest Bjgmf > 10 -17 — 10 -18 G in the local cascading 
environment. In contrast to these studies, we systemati- 
cally investigated effects of a wide range of uncertainties 
using detailed 3D Monte Carlo simulations to conclude 
that HO remains compatible with current observations. 

Two of these blazars RGB J0710+591 and 1ES 
1218+304, had VHE observations conducted during the 
Fermi mission and the measured VHE fluxes of these 
sources were assumed to be representative of the average 
VHE activity over the HE measurement period. For four 
other sources 1ES 0347-121, 1ES 1101-232, H 2356-309, 
and RGB J0152+017, the VHE component was measured 
prior to the Fermi mission and it was assumed to be rep- 
resentative of the average VHE activity during the Fermi 
observations. It was found that, for all of these sources, 
with Fermi-LAT Pass 7 data and a more general set of 
intrinsic source models of gamma-ray emission, HO can- 
not be rejected at > 95% confidence level f e yen w ith a 
typical EBL model, e.g. Domm guez et al.l (j2011f ) and 
models referenced within). 

A single source, the extreme blazar 1ES 0229+200, did 
allow for a rejection of HO at more than 99% confidence 
level, u nde r the standard assumptions as discussed in 
Section l4~3l However, if the EBL model SED is decreased 
in the visible and near infrared band to within the un- 
certainty limits of measurements based on the galaxy 
counts, the Bigmf lower limit cannot be maintained. 
HO can also be reconciled with the data if the VHE mea- 
surements of its flux conducted over the period of forty 
hours are not assumed to be representative of the past 
3.25 years of activity of this source. We must also point 
out that a thorough relative energy calibration of the 
Fermi-LAT and IACTs has not yet been reported and an 
error in the energy scale of VHE instruments may have 
significant impact on the conclusions regarding Bigmf- 

It should also be noted that a claim of the rejection 
of HO based on the observational data of a single source 
should not be considered robust due to possible varia- 
tions of the magnetic field in the local environment of 



the source. Figure [TTl shows the number density of the 
secondary electrons from 1ES 0229+200 for a zero mag- 
netic field as a function of distance from the source for 
electrons with energies > 100 GeV and > 1 TeV. It 
can be seen that a significant fraction of the secondary 
gamma-ray flux can be isotropized if cascading begins in 
the vicinity of several tens of Mpc of the source in the 
magnetic field of filaments which is expected to be or- 
ders of magnitude lar ger than the mag netic fi eld in the 
voids. Recent work bv lSutter et al.l (|2012f ) and lPan et al.l 
(2012) on identifying and characterizing the Voids in the 
Sloan Digital Sky Survey Data release 7 (|Abazajian et al.l 
(2009)), suggests that the volume occupied by them ac- 
counts for around 40 - 60 % of the total volume of the 
universe. It is therefore critical for studies of the Bigmf 
to know the exact line of sight distribution of voids par- 
ticularly in the few hundred Mpc vicinity of the source. 

Another avenue to reduce the secondary gamma-ray 
flux and reconcile it with observations and HO has re- 
cently been proposed by llBroderick et al.l (|2Q12h . The 
authors propose that beam-plasma instabilities could de- 
velop between the electrons and positrons in the cas- 
cade and electrons of the plasma in the voids. If such 
a collective interaction between two populations of elec- 
trons indeed exists, the energy dissipation rate of elec- 
tron positron pairs into modes of the plasma waves in 
voids may become significantly larger than the energy 
loss through IC scattering. Parameters of the plasma in 
the voids are largely unknown, but we will assume that 
the mass of the baryonic matter in the voids does not 
exceed the difference between the baryonic mass identi- 
fied through analysis of CMB fluctuations and the mass 
found in the galaxies and galaxy clusters. The cosmo- 
logical density of the universe is 5.9 baryons m -3 of 
which 4.6 % represents baryonic matter, which gives 
2.7 x 10 _T baryons cm -3 where 90% of these baryons 
were identified in filaments and clusters of galaxies fill- 
ing about half the volume of the universe. Therefore, 
the amount of baryonic matter in the voids cannot ex- 
ceed the remaining « 0.5 x 10 %, making the density esti- 
mate of the electron plasma in the voids n e « 1.4 x 10 -8 
cm -3 or lower. The plasma frequency of these elec- 

1/2 

trons, = 27r/p 5 e = (47rn e e 2 /ra e ) , is estimated to 

1/2 

be f PiG « 1 Hz (n e /10- 8 ) 1 . With no sources of en- 
ergy in the voids, such a plasma is thermalized with the 
CMB radiation on time scales of a few collision times, 
(^cmb^Tc) -1 , which is about 10 4 years. The screening 

1/2 

Debye radius of this plasma is = (&T/87m e e 2 ) = 

8.1 x 10 4 cmy/{T/2.7K) (10- 8 cm- 3 /n e ). We note that 
this screening radius is much less than the typical dis- 
tance between electrons in the beam, n^ 1 ^ 3 as estimated 
in our simulations shown in Figure [17j which is 10 6 — 10 T 
cm ^> A£>. It appears that electron positron pairs of the 
beam are nearly completely screened by the electrons of 
the plasma in the voids, and therefore, cannot interact 
collectively to generate a higher rate of energy transfer 
to plasma oscillations. 

Furthermore, one may argue that this screening can 
be avoided if the density of electrons in the voids is at 
least 2-3 orders of magnitude lower than the given upper 
bound estimate. For this regime, we consider a hydro- 
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Figure 16. RGB J0152+017 (a) (left) confidence level for a given cm, and e c . (b) (right) Simulated dFED of the best fit model of a = 
1.9, e c = 3.16 TeV, showing both the prompt and secondary cascade contributions to the total dFED, along with the HE and VHE data. 
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Figure 17. Electron Density as a function of distance from source 
1ES 0229+200, assuming an normal EBL model (EBL model 1) 
and a spectral index in the VHE band (a) of 1.3 and cutoff energy 
(e c ) 3.16 TeV. The black (solid) line traces the density of electrons 
with energy > 100 GeV, and the red (dashed) line represents the 
density of electrons with energies > 1 TeV. 



dynamical approximation of the beam-plasma instability 
for the two populations of electrons with zero tempera- 
tures and one beam at relative velocity c. The dispersion 
relation for these waves is given by 
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where oo p ^ is the plasma frequency of the electrons and 
positrons in the beam, and where we introduced the di- 
mensionless variables x = Re{uj)/uo p ^ e , y = Im(oo)/oo Pye , 

z = k • c/oo Pje and k = u) Pl b/w Pje = y/nb/n e <C 1. The 
real and imaginary part of this equation can be sepa- 
rated, and the unstable (y ^ 0) solution should satisfy 
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where e = \J zjx — 1. The solution exists for k < 
e < ft 1/3 and for < y/x < 3 1/2 (ft/4) 2 / 3 . Since 
k <C 1, this would require that \z\ / \x\ ~ 1 and 
that \z\ < 1. The last condition can only be sat- 
isfied for plasma waves with wave vector precisely 
perpendicular to the beam velocity. The cosine of 
this angle must satisfy the following condition cos# 
< (XD/coOp^kXD)- 1 = (k B T/2m e c 2 )^ 2 (k\ D )- 1 = 
1.5 x 10- 5 (T/2.7K) 1 / 2 (A:A D )- 1 , in which k\ D > 1 to 
enable collective interaction of plasma particles without 
significant screening. It is evident that maintaining such 
a requirement for the instability to grow is impossible 
even in the Biqmf = regime, because QED inter- 
actions introduce perturbations in the velocity distribu- 
tion of electrons of the beam of 7 6 _1 order, which for the 
bulk electrons of interest with energies larger than 100 
GeV, is comparable to or larger than the limit on cos#. 
Thus astrophysical inhomogeneities in the angular dis- 
tribution of the beam electrons should rapidly suppress 
the development of such instabilities. A detailed account 
of the relaxation of beam plasma insta bilities in cosmic 
voids has recently been investigated by iMiniati fc Elvivl 
(2012) with the conclusion that the relavistic pair beams 
of blazars remain stable on timescales much longer than 
the characteristic IC cooling time of electrons, and col- 
lective plasma-beam interaction effects in the voids are 
negligible. 

An alternative approach to reconcile a lack of the sec- 
ondary HE radiation in the observational data of ex- 
treme blazars is to suggest that a significant part of the 
VHE emission originates from proton-initiated cascades, 
rather than from direct photon s, if CR protons we re also 
accelerated by the same source (jEssev et al.l ()2Qllh ). The 
universe is nearly transparent to 10 16 - 10 19 eV protons 
on spatial scales of hundreds of Mpc and therefore, elec- 
tromagnetic cascades can be initiated by these particles 
in a random location along the line of sight. Crossing 
small regions of intense magnetic fields such as those 
present in clusters and filaments represent difficulties for 
this mechanism since they rapidly destroy the correla- 
tions between the cosmic ray directions and the source. 
Moreover, it is likely a formidable task to devise a mech- 
anism by which cosmic rays have been accelerated to ul- 
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tra high energies in the source by intense magnetic fields 
satisfying the Hillas condition and then highly collimated 
along a particular direction, eventually entering regions 
of potentially extremely small magnetic fields of voids 
without being significantly disrupted by the intermedi- 
ate magnetic fields. If this mechanism takes place in 
nature, it has two distinct observational characteristics, 
namely that the VHE radiation produced should show 
little evidence for variability and correlations with other 
wavelength bands such as X-ray, etc. and most impor- 
tantly, should be accompanied by higher energy gamma- 
rays, with energies above 10s of TeV in the VHE spectra 
of extreme blazars with z > 0.1. No such observational 
evidence has been collected so far to necessitate consider- 
ations of this scenario, which also lacks an explanation of 
the origin of the collimation. This m echanism has been 
recently studied (jMurase et al.l (j2Q12h ) in the context of 
1ES 0229+200 with the conclusion that the detection of 
larger than 25 TeV photons would provide an indication 
of acceleration of ultra high energy cosmic rays in this 
source. 

In summary, we have considered the observational data 
of 7 extreme TeV blazars and performed detailed simu- 
lations of cascading in the intergalactic magnetic field 
of the cosmic voids. We find for all sources, except for 
1ES 0229+200, no evidence to claim a lower bound on 
the IGMF or exclusion of HO (Bigmf = hypothesis), 
due to existing systematic uncertainties in source mod- 
els. For an explanation of the 1ES 0229+200 data set, 
it may also be necessary to account for astrophysical un- 
certainties such as the energy density of the EBL or the 
magnetic fields in the local (few hundred Mpc) environ- 
ment. 

The simulations were performed on the UCLA hoff- 
man2 cluster as well as the Joint Fermilab - KICP Su- 
percomputing Cluster. V.V.V. gratefully acknowledges 
grants from UCLA and S.P.W. gratefully acknowledge 
grants from Fermilab, Kavli Institute for Cosmological 
Physics, and the University of Chicago. This research is 
supported by the U.S. National Science Foundation un- 
der grants no. (PHY-0969948, PHY-0422093, and PHY- 
0969529). 
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